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The reconstruction of the CMB lensing potential is based on a Taylor expansion of lensing effects 
which is known to have poor convergence properties. For lensing of temperature fluctuations, an 
understanding of the higher order terms in this expansion which is accurate enough for current 
experimental sensitivity levels has been developed in Hanson et. al. (2010), as well as a slightly 
modified Okamoto and Hu quadratic estimator which incorporates lensed rather than unlensed spec- 
tra into the estimator weights to mitigate the efi^ect of higher order terms. We extend these results 
in several ways: (1) We generalize this analysis to the full set of quadratic temperature/polarization 
lensing estimators, (2) We study the effect of higher order terms for more futuristic experimental 
noise levels, (3) We show that the ability of the modified quadratic estimator to mitigate the effect 
of higher order terms relies on a delicate cancellation which occurs only when the true lensed spectra 
are known. We investigate the sensitivity of this cancellation to uncertainties in or knowledge of 
these spectra. We find that higher order terms in the Taylor expansion can impact projected error 
bars at experimental sensitivities similar to those found in future ACTpol/SPTpol experiments. 



I. INTRODUCTION 

Over the past year, data from two ground based tele- 
scopes, ACT and SPT, have resulted in the first direct 
measurement of the weak lensing power spectrum solely 
from CMB measurements [1, 2]. In the coming years, 
the data from Planck and upcoming experiments ACT- 
pol and SPTpol will begin probing this lensing at much 
greater resolution. The state-of-the-art estimator of weak 
lensing, the quadratic estimator developed by Hu and 
Okomoto [3, 4], works in part through a delicate cance- 
lation of terms in a Taylor expansion of the lensing effect 
on the CMB. In this paper we present a simulation based 
approach for exploring the nature of this cancelation for 
both the CMB intensity and the polarization fields. In 
particular, we study a slightly modified quadratic esti- 
mator which: incorporates lensed rather than unlensed 
spectra into the estimator weights to mitigate the effect 
of higher order terms; and uses the observed lensed CMB 
fields to correct for the, so called, Ni^^ bias. 

The simulation methodology presented here allows a 
stochastic exploration of the higher order bias terms of 
the quadratic estimate and can be used to reduce the 
computational load associated with iterative de-biasing 
algorithms for the quadratic estimate. In this paper, 
we use our simulation methodology to present a de- 
tailed study of the, so called, Ni^^ and Ni"^"^ bias for 
the full set of quadratic temperature/polarization lens- 
ing estimators. The N^^'' bias was first explored for 
the standard flat sky quadratic estimate in Kesden et 
al. [5]. For full sky CMB temperature maps, Hanson 
et al. [6] developed an approximation to the higher or- 

der bias terms, including TVj ' , which is accurate enough 
for current experimental sensitivity levels as well as for 



* Supported by NSF grant 1007480.; anderes@stat.ucdavis.edu 



the slightly modified quadratic estimator which incorpo- 
rates lensed rather than unlensed spectra into the esti- 
mator weights. We generalize this analysis to the full 
set of modified quadratic temperature/polarization lens- 
ing estimators and demonstrate that, indeed, the lensed 
spectra weights mitigate the combined higher order bias. 
However, this mitigation is obtained only by an increase 
in the magnitude of both A^/^^ and Ni'^^ to the extent 
that they nearly cancel. We explore the extent with 
which this cancelation is sensitive to fiducial uncertainty 
in the way the lensed spectra weights are computed. 
We find that, under experimental conditions similar to 
those in future ACTpol and SPTpol experiments, the 
EB quadratic estimator is not sensitive to low I fiducial 
uncertainty whereas the EE and TE are sensitive to the 
point of degrading inferential power. 

The remainder of the paper is organized as follows. In 
sections II, III and IV we give an overview of the modi- 
fied quadratic estimate, derive the spectral density of the 
quadratic estimate in terms of higher order bias terms, 
and discuss the estimation of the spectral density of the 
lensing potential. In Section V we present two simula- 
tion based methods for estimating the higher order bias 
terms. The first simulation method works exclusively for 
estimating A^/^\ and is mainly used to validate the sec- 
ond algorithm which can produce all higher order terms 
N^^^ for j > 1. In Section VI we use these methods to 
study the higher order terms for experimental noise lev- 
els similar to those found in future ACTpol and SPTpol 
experiments. The paper concludes with the Appendix 
which gives fast Fourier transform (FFT) algorithms for 
computing the modified quadratic estimate. These algo- 
rithms extend the FFT techniques developed in [3] to the 
computation of all quadratic normalization constants and 
provides fast (non-stochastic) algorithms which extend 
the simulation techniques found in [1] for computing the 
'Gaussian bias' from the lensed CMB four-point function 
for all temperature/polarization quadratic estimators. 
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II. THE QUADRATIC ESTIMATOR 

The effect of weak lensing is to simply remap the CMB 

temperature T{x) and Stokes polarization fields Q{x) 
and U{x) for a flat sky coordinate system a; G M^. Up 
to leading order, the remapping displacements are given 
by V(/)(a;), where (f>{x) denotes a lensing potential and 
is the planar projection of a three dimensional gravita- 
tional potential (see [7]). Therefore, for any CMB field 
X € {T,Q,U} the corresponding lensed field can be writ- 
ten X{x + V4>{x)). For the remainder of the paper we 
let 

X{x) = X{x + V(l){x)) + N'^ix) 

denote the corresponding lensed CMB field with addi- 
tive independent experimental noise given by (which 
includes a beam deconvolution) . Using this notation 
the corresponding lensed E and B modes are given by 
El = - coa{2Lpi)Qi - sm(2(pi)Ui and Bi = sm{2ipi)Qi - 
cos{2ipi)Ui where Xi = J ^e^^^'^X{x) (unitary angu- 
lar frequency) and ipi denotes the phase angle of fre- 
quency I. 

For any field X, F e {T, Q, [/, E, B) the spectral den- 
sity Cj^ is defined to satisfy {XjYi,) = Si^iiCf^ where 



The angled brackets (•) denote ensem- 



ble averaging (or expected value) over both the CMB 
fields and the large scale structure given by (f). In ad- 
dition, we let {•) xY denote expected value with respect 
to the unlensed CMB fields T, Q, [/, E, B and {^^ denote 
expected value with respect to large scale structure given 
by (j). Throughout this paper we stipulate (j) is indepen- 
dent of T,Q,U,E,B which implies: (•) = {{■)xy)<p- We 
let denote the lensed CMB spectral density without 
experimental noise and let C/^^ = Cf^ + Cj^ ^ and 

Cj-^p = + Ci^^^^ denote the corresponding un- 
lensed and lensed spectral densities with the additional 
experimental noise. 

The quadratic estimate, based on two lensed CMB 
fields X and Y, is derived from the following two state- 
ments: 

(Xk+in)xr'^^lfl%' whenZ^O; (1) 
(X^+in) =SiCiZ.^ (2) 



which hold for any X,Y € {T,E,B} and where the co- 
efficients are given in the Appendix. Equation (1) 
approximates the cross frequency correlation (at separa- 
tion lag I) induced by the nonstationaritj^ in X (when 
regarding 4) as a fixed nonrandom field). This is derived 
through a Taylor expansion of the lensing operation for 
any X G {T, Q, U} 



X{x) = 5^X{x) + S^X{x) + S'^X{x) + ■ • 



(3) 



where S°X{x) = X{x) + N'^ix), S^X{x) = 
VX{x)'Va(l){x), etc. When X e {E,B} one de- 
fines 6^ El = -cos{2(pi)6^Qi - shi{2ipi)5Wi and 5^Bi = 



sin{2(pi)5^ Qi—cos{2ipi)S^Ui. Then by expanding X^_^iYi^ 
with (3), regrouping terms by the order of (p, one obtains 
fi^cl^i = {5'Xk+i6°Y^)xy + {S°Xk+i6^Y^)xv which gives 
approximation (1). Equation (2), on the other hand, is 
obtained by treating both the CMB and the large scale 
structure cj) as random so that X, from this viewpoint, is 
isotropic (but non- Gaussian). 

Hu and Okamoto [3, 4] used approximations (1) and (2) 
to construct the optimal quadratic estimate of cj) based 
on X and Y as follows 
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(4) 



where = 2nf{^[C^i^^^^C^^J \ The normalizing 
constant Af^ is determined through an unbiased con- 
straint. In particular, using the fact that //^ is real we 

have that = ^i^r I j^mi'^r+u^cz.. 

by equation (1). Then requiring that (0f^) = <pi deter- 
mines Af^ as follows 



J w'^' 



k\ ^k+l,e^p^k,e^p 



(5) 



A. Lensed versus unlensed weights 



There is a small modification to the standard quadratic 
estimate (pf^ which can mitigate the low I bias (arising 

(2) 

from the A^,^ ^ t( Tm discussed in the next section) when 
using the observed power \4>f^\'^ to estimate Cj*"*". This 
modified estimate, denoted , is obtained by replacing 
all occurrences of unlensed spectra in gf^ and Af^ with 

the corresponding lensed spectra. In particular ^ 
defined as 



IS 



^Af 



27r 



(6) 



XY 
I ' 



Where Af- = [j ^\~gi%\'C^+i,^^C^^J and ^i^l is 
obtained from gf^ by replacing every occurrence of C, 

C/^^ and C^^ with the corresponding lensed spectra 
C^-^ and Cp' . For example, when X = Y — T one has 

^^_i-[{k+i)ciii-kcr] 

9i,k = ' 



'-'fc-l-I.oxp^fe.cx 



9i,k 



l-[{k + l)C, 



TT 

k+l 



^k+l,e^V^k,e^V 



Notice that the estimate (pl^ is not normalized to be 
unbiased. Indeed from equation (1) one has 



d?k 



(2^) 



~XY „XY /--»XX /^Y 

2 9l,k9l,k^k+l,e^p^k 



YY 

exp 



3 



III. THE SPECTRAL DENSITY OF THE 
QUADRATIC ESTIMATE 

In this section wc dciiivc! the fohowing all order decom- 
position of the spectral density of the quadratic estimate 

(0r <^r* ) = Si-i' [Cf^ + + nI'^ +■■■] (7) 

which will then be used, in a subsequent section, to derive 
the estimation bias for Cf"^. The first term A^/^-* is related 
to the disconnected terms of the lensed CMB four-point 
function, whereas the higher order terms A^/"'^ for j > 
1 arc related to the connected terms of the four-point 
function segmented by the order of Cf"^. Most of this 
section focuses on the quadratic estimate followed 
by a brief discussion of the corresponding decomposition 
for the modified quadratic estimate (f>^ . 

Our derivation of the spectral density of in the 
flat sky is similar to the analysis of the full sky trispec- 
trum done in Hanson et al. [6]. One starts by relating 
((^^0^') to the lensed CMB four-point function by dis- 
tributing the expected value as follows 

X /^/ ^9u.9p:k'{X,+iYX+i'yk'). (8) 

To decompose (8) one then expands the four-point prod- 
uct term in the above integrand by expanding the lensed 
CMB Taylor expansion (3) to obtain 

^ ^ ' 

decomposes further into connected 
and disconnected terms 

(9) 



A. Disconnected terms 

After distributing the expected value in the right hand 
side of (9) through the j-fold convolution which makes 
up 6^Xi and subsequently applying Wicks theorem, 
one can further decompose 1^5' X,^^i6^Y^5PX^,^i,5'^Yk') 
into what are called connected and disconnected terms. 
The disconnected terms in the four-point product are 
the terms which factor into cross-spectra of the fields 
5'X, S^X, SPX and S'^X. For example, if {i,j,p,q) = 
(1,1,0,0) and (p is assumed independent of X and Y 



then 

{5'X^^i6'Yj:S°Xl,^i,S°Yk>) 

I 1 I 1 

= (^S^Xi^_^i5^Y^ (5*Vfc/ ) -f- disconnected term 



+ (^5^X f^_^i6^Y^ S^X^,_^i, S^ki ) <— connected term 



+ (^6^X i^_^i5^Y^ 5'^X^,_^i, ^'Vfe/ ^ 4— connected term 



where the top contraction symbols correspond to pairing 
the CMB fields and the bottom contraction symbols cor- 
respond to pairing the lensing potential in 6^X and S^Y. 
The contraction pairings on the above disconnected term 
results in a product of two spectra as follows 

{6'x'^Y:6'^xi,^,S%.) 

= {S'X^+iS'y:){S'XI,^, 6"V^,}. 

In a similar manner, for general {i,j,p,q), the discon- 
nected terms can be grouped into the three types: one 
for each possible configuration of the of top contraction 
symbols. Then regrouping all disconnected terms in (9), 
by top contraction type, results in the following three 
terms: 

disconnected terms in (^S^X|^_^_^6■'Y^6PX^,_^_l,S'^Yk') 



— {^k+iYk ^k'+i' Yk' ) + {Xk+iYk ^k'+i' Yk' ) 

I 1 I 1 

+ {^k+iYk ^k'+i' Yk' ) 

^ V ' 

=0 when I or I' 

= (^k+l,e.pCk^e.Jk+k'+l'Sk+k'+l (10) 
+ (-'k+l,e^p'^k,<,^p(>k-k'dk+l-k'-l' 

where the last line is obtained by applying approximation 
(2). Substituting the four-point product term in (8) with 

(10) results in, what is typically called, the 7Vj^°^ bias: 
disconnected terms in {4>i^4>P^ ) = ^l-l'^i^^ 

where 

Ar(0)_rjxyi2 f d k f „XY P^XY P^XY 

^^l = \ J V '^fe-|-J,exp<-^fe,exp 

I I XY \ 2/^lXX p^iYY \ 

'T~\9l,k\ '-"fc-|-J,exp'-"fc,exp J ■ 

B. Connected terms 

The connected terms decompose further into what wc 
call the 'first connected terms' and the 'second connected 
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EE estimator 

— l^Ct'/i 

estimated l^Cf"^/4 (lensed weights) 
estimated l^Cf'^/i (unlensed weiglits) 
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EB estimator 






* estimated l'^Cf*/i (lensed weights) 

O estimated (unlensed weights) 

* '*(iV,''jj, )/4 (lensed weights) 
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FIG. 1. These figures demonstrate the difi'erence between EE and EB quadratic estimators when using lensed versus unlensed 
spectra into the estimator weights. The main feature of the left hand plot is the reduction of bias — comparing with 'o' — in 
the EE estimator when the lensed spectra are used for the estimator weights. Also by comparing with 'o' and with 
'o' it is clear that this reduction is obtained by an increase in magnitude of both the so called A'^'^' and A'^' bias to extent 
that they nearly cancel. The first order bias term l'^{Ni]^^^) /A is denoted for (^f^ and 'o' for ■ The second order 
bias term ^^(-Aj''^'j^^)/4 is denoted for (f>f^ and 'o' for 4>i^ ■ Finally, the expected value of the spectral density estimates 
l^{So^\4>r'\^ - KoL>/4 and l*{So^\4>r'\^ - A^/°oL>/4 are denoted by 'o' and respectively. See Section VI for details. 



terms'. There are only four 'first connected terms' and 
are defined as follows: 

first connected terms m ^ ^ (^S^ X j^^iS-' X^, 6''Yf,'') 

^ {S' Xk+iS''Y^6'Xp ^i,S^Yk,) 

+ {d' Xk+iS''Y:s"Xl,+i,S'Y ^,) 
+ (S^X^Y^S^xf^Y^,) 
+ {d'Xk+iS' Y^5°xi+i,S'Y ^,) 

Then, by substituting the four-point product term in (8) 
with the first connected terms, one gets 

first connected terms in (j)^ ^ = Si^irCf'^ . (11) 

The remaining connected terms in the four-point prod- 
uct of (8), called 'second connected terms', are then re- 
grouping corresponding to the order of cf). After noticing 
that any term of order 0^-'+^ has expected value zero one 
obtains the following expansion 

second connected terms in(^<l)f^ (p^ ^ 

where A^j'-"''' is of order <jp^ . Putting all disconnected and 
connected terms together gives the desired expansion (7). 



For reasons which will become clear in the next section, 
we define A^/^^ in a slightly different, but equivalent, way 
that extends more naturally to the modified quadratic 
estimate. In particular, instead of defining A^;*'"'^^ as the 
as the total contribution of the second connected terms 
in (8) of order (/)^, we define Ni^^ simply as the total 
contribution of all connected terms of order (ji^ minus 

C. The lensed weights 

For the modified quadratic estimator, 4)f^ , one can 
derive the expansion (7) with a few minor adjustments. 
In particular, the disconnected terms can be written 

disconnected terms in (^4>i^ (j>v^ ) — Si-i'N^^"^ 

where 

Ar(0) _ r 4X1-12 f d k / 

,cxp fc,cxp 

+ 1 9iM I ^ C'fe+i , cxp C'fe^xp ) • 

The first connected terms are slightly different due to 
the intentional bias in the modified quadratic estimate 
yielding 

first connected terms in{4>^4>^') = 5i_i,[B^fCf''' . 
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FIG. 2. The first order bias term ^''(A^/,;i,)/4 is denoted for (^f^ and 'o' for 4>f^ . The second order bias term l*{Nj; J^^)/A 
is denoted for and 'o' for cf)f^ . Finally, the expected value of the spectral density estimates /*{(5(^^|(^ 
and Z*{(5(7^|0f^|^ — -^/oL)/'^ denoted by 'o' and respectively. See Section VI for details. 



xy\2 
I 



r-<l>/4 



After these adjustments are made the remaining terms 

Ni are defined exactly the same way as for (f>P': when 
(i) 

j > 2, Ni is defined as the total contribution of the 

connected terms in (8) which are order 0-'; and N^^^ as 
the total contribution of the connected terms in (8) which 
are order tp"^ minus Cf"^. This is consistent with the pre- 

vious definition of ' for and preserves the natural 
expansion 

(^r ) = 5i-v [Cf^ + + iv/'^ + • • • ]. (12) 

IV. ESTIMATION OF Cf 

In this section we show how the expansions (7) and (12) 
are used to derive and analyze estimates of Cf"^. All the 
estimates presented in this section can also be addition- 
ally radially averaged, with inverse variance weights, to 
reduce estimation variability. We derive the following re- 
sults using the estimate (j)f^ and simply remark that the 
results can be similarly derived for the modified quadratic 
estimate (j)f^ . 

For current experimental conditions the first term N^'^'' 
dominates the sum J^'jLq ■ Therefore a natural bias 
corrected estimate of Cf'^ is given by SQ^\4>f^\'^ — A"/""*- 
Notice, however, that to compute N^^"^ one needs a model 
for the lensed spectra with experimental noise: C/^^^, 
Cr^ and CP' . This normally requires knowledge of 

t , GXp t ,CXp i/ i O 

the very quantity we are estimating: Cf"^. One way to 
circumvent this difficulty is to replace the experimental 
lensed spectrums with estimates from the observations 



X and Y . This results in the observed Aj bias and is 
given by 

+ \ga\'\^k+i\'\n\'). 

Using nI;11^ in place of Af/°^ yields S^^\(j)f^\^ - A^/°|,^ as 

an estimate of which docs require knowledge of Cf"^ 
to compute it. 

Remark: Up to this point we have been assuming 
an infinite sky when computing the Fourier transform 
f However, the above Fourier transforms will 

typically be done on a pixelized periodic finite sky. In this 
case, Sq is approximated as 1/AZ where Al is the area el- 
ement of the grid in Fourier space induced by finite area 
sky. For the remainder of the paper we do not distin- 
guish the finite versus infinite case and simply equate Sq 
with l/Al leaving it understood that equality holds in 
the limit as AZ — > 0. 

A. The bias of S-^\4>n'' - <oL 

We will derive three main expansions, which decom- 
pose the bias: {S^^l^ri^ " ^i^L) " C*/"^ (and similar 
definitions for the modified quadratic estimate). These 
expansions will be denoted as follows 

(^f = ^/"^ + + ^/'^ + • • • (13) 

\4,f^\^ = ef^ + e^p + + . . . (14) 
<L = ^^''^H-^'^^ + -4^'' + --- (15) 
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'''(-'Vi't'l,. >/* (lensed weights) 
'■''('^1%, >/4 (unlensed weights) 
1*{nI^^^^)/4 (lensed weights) 

>/4 (unlensed weights) 




600 800 
wavcnumb cr 



FIG. 3. The first order bias term /*(iV/,tL)/4 is denoted 



for </>f and 'o' for 



, aiiu ^ ioi (/)f^. Tlie second order bias term 
l'^{N^%^)/4 is denoted for and 'o' for 0f^. Fi- 

nally, the expected value of the spectral density estimates 

iHSo'lin' - <L>/4 and l^So'l^n' - <L)/4 are de- 
noted by 'o' and respectively. See Section VI for details. 



we have that 



(0)N 



(18) 



where N^^^^ is of order c/P^ . Therefore the order (j)'^^ term 

in the bias of S^^l^f^]^ - 7V/°|, is given by (N^'I ) for 
J >1. 

Remark: Notice that the zero order term can be com- 
puted easily 



(0)N 



(2^)2 

^i/2,oxp^i/2,oxp + l'-'!/2,cxp7 



when 1^0. Moreover, since this term does not depend 
on (j) one can simply subtract it out of the estimator when 
using (5(7^|<^f^P — -^/oIb for estimation. 



B. The difference between Nj; ' and (A^/ Jb^) 



where the terms G'l''^ and are each of order 

j in (j). Expansions (13), (14) and (15) are all obtained 
by expanding and Yf, using (3) then regrouping the 
terms by the order of (/>. This expansion and subsequent 
re-grouping yields the following analytic expressions for 
each term; 

n ^ 



p=0 
j 



p=0 



By matching the right hand side of equation (14) 
with the right hand side of (7) one gets that SqN^^^ 



equals 



?(2)\ 



minus SoC[ and any disconnected terms. 



In contrast, <5o(A'j'^|^^) equals 



^ SqCi and 



{Sq^i'"^^). Therefore the difference between ^o-^/^'' and 
'^o(-^/oL) is difference between (i5o=y^^^^) and the 



(2) 

disconnected terms in (^j ). By expanding SqC^^^ 



J2'^=o{^i^ definition of SgN^^^ (which equals 
all the disconnected terms), and retaining only the order 
(j)^ terms (what remains equals all the disconnected terms 

in (^f ^)) we get 



1 9i%9i^k-i 



^xv 
k+l 



^XY,2-p 



\9i,k\ '^k+i ^k 



where ^^^^ = ELoS^XkS^-PY^ so that XkY^ = 



EOO 
p^O ^ k 



Using the fact that 
one then gets that 



+ \9rk\\^kin{^v^'-'))- (19) 

Therefore the difference between SoNj^^^ and ^o(-^/oL) 

given by the difference between (19) and ((^o^ ) which 
equals 



(2i) 



Now by defining A'^^^^j^^ for j > 1 as 



N, 



Z.obs 



-5o^(^^")). (16) 



when j = 1 
when j > 1 



(17) 



5^.P tf0^'2-j 
k 



Y^yiArrJ ^,{9^k9r-k-i {^^^r^, 

+ \9aH^k+f^r''-'))- (20) 



p=0 



For the experimental conditions analyzed in this paper 
this difference is small. 
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V. FAST MONTE CARLO ALGORITHMS 

In this section we give two simulation based methods 
for quickly estimating N^-'^ and {Ni''^^^^) for j > 0. The 
first method simply observes that each term ff'i ^ 

and ■-yVi'"''^ have fast Fourier transform characterizations 
which can be nacA for simulating N^''^^^^ and by averag- 
ing multiple realizations — for estimating (-/v/^^bs)- This 

(7) 

algorithm also extends to A^j by replacing the expan- 
sion of Ni^]^^ in (15) with the corresponding expansion for 

Ni^\ The second method is exclusive to Ni^'' and uses 
correlated and uncorrelated CMB fields to mimic the ap- 
propriate Wick contractions for the connected terms in 
(8). 



(5q "'^l^f^P — -^/^L ^ where /ij is an approximate bias 
correction based on the estimates of the higher order bias 
correction terms {N^^^^^y Notice that for fast re- fitting 
algorithms, it maybe be advantageous initially to toler- 
ate relatively large the Monte Carlo error bars, then it- 
eratively increase the number of Monte Carlo samples at 
each re- fit. 

We finally mention that the above simulation methods 
which are used to estimate (-/v/'^^^) can also be used to 
estimate Np\ The only change is to replace each term 



k+l ' ' 



^k+f ^^'^ S^^''^^ ^ in the definition 



of .y^i'^^^ with their respective expected values: {^^i), 

( <^-,2,-p^^ (^-f ) and 



B. Coupling lensing fields for iVj 



(1) 



A. FFT algorithms for Nl 



U) 



Ob. and iV, 



0) 



The fast simulation techniques presented in this section 
depend on the fact that the transforms which character- 
ize (j)^ and (f>f^ can be derived as Fourier transforms of 
point-wise products of gradients in pixel space. This was 
first utilized in [3] and [4] for the flat sky quadratic es- 
timators. In the appendix, we present these transforms 
along with some additional FFT transforms which allow 
fast simulation of the fields ^l^^ and This is 

the basis of the algorithm which then uses equation (17) 

to simulate N, ' . 

For ^/■'^ and ^j'''' first notice that each term 6^X1 can 
be easily simulated in the pixel domain since S^X{x) is 
the point-wise product of derivatives of X{x) and (j){x). 
Moreover the quadratic estimate applied to these terms, 
resulting in / ^5/^ [5PXk+iS^-PY^] , is also easily 
computed by direct application of the formulas presented 
in the Appendix. Now summing over p G {0, 
gives fast simulation of ^/''^ which, in turn, gives ^i''^ 



(0) 

I ' • 



Fi- 



by taking quadratic combinations of 

nally, to simulate -y^i^^^ start by noticing that each term 

^k^'^ = Ep=o SPXkS^-Py* simulates easily from SPXk 

and 6^~PY^. Then to recover c/^^"*^ one uses the FFT 
transformations presented in the appendix for quick sim- 

ulations of [Ar? l£^9rk9r-k-i ^^T+l ^r'" and 

[4x^12 r_dik_\„xY\2 agxx,p ^YY,2-V 
L^i J J (27r)2 li/i,fe I '^k+l ^k 

Monte Carlo averaging n independent simulations of 
-^/IL ^^^^ yield estimates of (-/Vj^^^^) with error bars that 
can be approximated by si / ^/n where si denotes the stan- 
dard deviation of the samples at each frequency I. In ad- 
dition, when the noise and beam structure are isotropic 
one can reduce the error bars by radially averaging each 
simulation of N, , . These error bars can then be used 

t ,ob.s 

for re-fitting algorithms where Cf'^ is iteratively fit to 



A second algorithm, for Monte Carlo estimation of 
Ni^\ is to use correlated and uncorrelated CMB fields 
to mimic the Wick contraction structure appearing in 
the definition of Ni^\ The algorithm is easiest to illus- 
trate with the EB quadratic estimator since the defini- 
tion of Ni^^ only involves a small number of connected 
terms in (8) (when assuming a zero B mode in the un- 
lensed CMB polarization field). The algorithm is derived 
by first noticing that the sum of the connected terms of 



order 



in (0f ^ 



) is given by 



Ol-l' 



[Nr + en 



(Pk 



I — 



XY XY 

9i,k9i'k' 



{5^Ek+i5^Bl5''El,^i,5^By) 



These two Wick contraction terms have a Monte Carlo 
characterization as follows. Let (Qj, C/j) and (Qj, U[) de- 
note two independent realizations of the CMB polariza- 
tion. Now let 0™ and (/)f| denote the EB quadratic 
estimator applied to the pairs {Ei,d^Bi) and (i?j, (5^i?j), 
respectively, where 5^B and 6^B' use the same lensing 
potential (/> so that 

6^Bi = smi2^i)S^Qi - cos{2^i)5^Ui 
S^B'i = sm{2ipi)6^Q[ - cos{2ipi)5^U[ 

where S^Q {x) = V''Q{x)V a<l>{x) and 5Ki'{x) = 
W Q' {x)V a4'{x) (similar definitions for 5HJ {x) and 
5HJ'{x)). Notice the independence structure of the sim- 
ulated fields-that E and B are independent of E' and 
B' — implies 



C2) 



ArAr 



27r 



XY XY 

yi,k Sl'.k' 



X {6''Ek+i6^Bl5°El,^i,S'Bk'). 
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FIG. 4. Spectral estimation bias averaged over the multiple bin I £ [10, 200] as a function of fiducial uncertainty. See Section 
VI for details. 



Similarly let <?!)f 3 and (j)f1 denote the EB quadratic esti- 
mator applied to the pairs and re- 
spectively (again using the same lensing potential 
Then 

2-K I 2n ^'''»^''>'»' 

I ^ I 

This leads to the following Monte Carlo averaging char- 
acterization of Ni^^ 

N^'^ = So'{4>ri4>f$' + 4>fX' ) - C't'^- (21) 

The above formula also holds for the lensed quadratic 

estimator t/i™ as well. It is easy to sec that this method 
can be extended to all other polarization quadratic esti- 
mators, but with decidedly more connected terms. 

Notice that there is an additional simplification when 
using the quadratic estimator, (pf^, without lensed 
weights. In particular, {^f^^ff ) = Cf"^ so that 
equation (21) simplifies to Nj;^^ = So\4>f^j^4>f^' ) . How- 
ever, this simplification docs not hold for (j)f^ since 
the bias factor Bf^ , defined in Section II A, implies 

One advantage of using this coupling technique is that 
the simulation of the and d'^Y can be done without 
including the additive experimental noise. To see why, 
notice that any CMB Wick contraction connecting 
and (which depends on the additive experimental 
noise) must have bottom contraction symbols connect 
S^X and 6^Y. This must yield a disconnected term which 
does not contribute to xj^K 



VI. SIMULATION 

We perform two simulation experiments under experi- 
mental conditions similar to those found in future ACT- 

pol/SPTpol experiments. The first simulation explores 

the bias terms X^-''' and (A^/^j^^) for the two quadratic 

estimators and 0f^. The results are summarized in 
figures 1, 2 and 3. For these simulations, to compute the 
lensed weights in (j)f^ , we use the same fiducial model 
for Cj^* as the simulation model for C/"^. In contrast, 
for the second set of simulations we explore the effect of 
uncertainty in the fiducial model for Cf"^ when comput- 
ing (pf^ . The results are summarized in figures 4, 5 and 
6. The main conclusion of the first set of simulations 
is that although (j)f^ does reduce low I estimation bias, 

this is accomplished by increasing X^^^ (or {Xi]^^^^') in 

the case that one uses 7V/"|,^ instead of xl;°^) to the point 

of canceling with Aj'^^ (or {X^'^J^^^') as the case may be) 

when the correct model for Cf"^ is used to generate the 
lensed weights. The second set of simulations show that 
this cancelation can be sensitive to the fiducial model for 
Cf"^ depending on which estimator one uses: TE and 
EE arc most sensitive, EB is least sensitive. At the end 
of this section we discuss the inferential implications for 
future experiments. 

The cosmology used in our simulations are based 
on a flat, power law ACDM cosmological model, 
with baryon density fib = 0.044; cold dark matter 
density ficdm = 0.21: cosmological constant density 
flA = 0.74; Hubble parameter h = 0.71 in units of 
lOOkms"^ Mpc~^; primordial scalar fluctuation ampli- 
tude Asik = 0.002 Mpc^i) = 2.45 x 10"^; scalar spectral 
index ns{k = 0.002 Mpc^^) = 0.96; primordial helium 
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^^Q-a TT estimator bias relative to fiducial Cf error ^^^-e TE estimator bias relative to fiducial error 




-0.3 -0.2 -0.1 0.1 0.2 0.3 -0.3 -0.2 -0.1 0.1 0.2 0.3 

(fiducial Cf^-true C**)/ (fiducial Cf*) (fiducial Cf^-true C**)/ (fiducial Cf*) 



FIG. 5. Spectral estimation bias averaged over the multiple bin I £ [10, 200] as a function of fiducial uncertainty. See Section 
VI for details. 



abundance Yp = 0.24; and rcionization optical depth 
Tr = 0.088. The CAMB code is used to generate the 
theoretical power spectra [8]. 

To construct the lensed CMB simulation we first gener- 
ate a high resolution simulation of 6 (a;) and the gravita- 
tional potential (l){x) on a periodic 25.6° x 25.6° patch of 
the flat sky. The lensing operation is performed by taking 
the numerical gradient of (f>, then using linear interpola- 
tion to obtain the lensed field &{x + Vc6(a;)). We down- 
sample the lensed field to obtain the desired 1.5 arcmin 
pixel resolution for the simulation output. Finally, the 
observational noise (with a standard deviation of 5/iK- 
arcmin on T and y/2 x 5/LtK-arcmin on E, B and Gaus- 
sian beam FWHM=1.5 arcmin deconvolution) is added 
in Fourier space. For all of the simulations we assume 
a zero B mode and a lensing potential </> which is un- 
correlated with the CMB. In contrast to the full lensing 
simulation, the pertabivc expansions given in Section V A 
only require simulation of unlensed CMB fields at the low 
resolution 1.5 arcmin pixels. 



Figures 1, 2 and 3 

Each plot in figures 1, 2 and 3 correspond to a different 
quadratic pairing X,Y & {T, E, B} and shows the Monte 

Carlo approximations to Z^(A^/_tL)/4 and ^^(A^/.!L)/4 
along with the expected value of the spectral density 
estimates n5^'\^r'? " <L)/4 and ;4(<5-i|<^-|2 _ 

-^/oL)/'^ over different realizations of (j), the CMB and 
the observational noise. Although not shown, the bias 
terms l^N^^'^ /A were also computed, using the coupling 
technique given in Section VB, and resulted in very simi- 
lar plots (mostly indistinguishable above the Monte Carlo 
error). The spectral density estimates are computed from 



the all-order lensed simulations whereas the bias terms 
are computed using perturbative expansions discussed in 
Section V A. The Monte Carlo approximations are based 
on 2000 independent realizations for the TT, EE and EB 
estimators and 18000 independent realizations for the TE 
and TB estimators. These estimates are then radially 
averaged on sliding concentric annuli with wavenumber 
bins of width 20 to yield the plots shown in figures 1, 2 
and 3. 

The main feature in these simulations is the large in- 
crease in (A'j' bias at low / for <j)f^ as compared to the 
corresponding quantity for (pf^ , especially for the EE, TE 

and TT estimators. In contrast, {N^"^^^^ also increases in 
magnitude but to a lesser extent, enabling the cancela- 
tion with (A^/^ti.)- Since the terms (A^/.^L) and 
are not individually small but instead cancel, there is the 
potential for this cancelation to be offset when there is 
uncertainty in the fiducial model for Cf"^ used to compute 
the lensing weights for the estimate ^f^. In the next sec- 
tion we explore this sensitivity by analyzing the resulting 
estimation bias as a function of fiducial sensitivity. 



Figures 4, 5 and 6 

To explore the effect of fiducial uncertainty when com- 
puting the lensed weights in (pf^ we fixed a fiducial model 
Cf"^ used to compute the lensed weights, then analyzed 
simulations under perturbations of Cf'^. In particular, 
we considered simulation models which differ from the 
fiducial model by a maximum of 35% only in the mul- 
tipole range [10,200]. The simulation models are of the 
form TiCf''' where Tj = 1 when / ^ [10, 200] and Tj = c 
when I S [10,200] where c ranges from 1.35 to 0.65. For 



each scalar c wc simulated 200 different CMB, lensing 
and noise fields. At each simulation we recorded the es- 
timation error, namely l^{6o^\4>r'? - -^£L " TiCf"^) /4, 
averaged over all frequencies in the / bin [10,200]. This 
error is then averaged all 200 simulations to estimate the 
bias. This bias is then plotted for each quadratic pairing 
X,Y £ {T, E, B} in figures 4, 5 and 6. The estimation 
bias is shown as '— •— ' with la monte carlo error bars 
attached as a function of the fiducial uncertainty 1 — c in 
the bin / G [10, 200]. The dashed line shows the standard 

deviation of l'^{5o^\(j>f'\^ - iVi°ij/4 and the solid black 

line plots the fiducial error (l-Tj)/4cf'^/4 averaged over 
I e [10,200]. 
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bias can nearly ignored completely. 



<PI 


~ bias range 


~ full sky 4a 


fiducial error (15%) 


TE 


±0.70 


±0.76 


±2.02 


EE 


±0.53 


±0.62 


±2.02 


TT 


±0.16 


±0.86 


±2.02 


TB 


±0.10 


±2.33 


±2.02 


EB 


±0.01 


±0.67 


±2.02 



TABLE 1. In units of spectral power, per 10~ , the above table 

shows bias range, approximate full sky standard deviation and 
fiducial error (at 15%) for each quadratic estimate (^™, i^f^, 
4>T^, 4>™ and 0f^. 



^^qS TB estimator bias relative to fiducial error 




-0.3 -0.2 -0.1 0.1 0.2 0.3 



(fiducial Cf*- true Cf *)/(fiducial Cf*) 

FIG. 6. Spectral estimation bias averaged over the multiple 
bin / £ [10, 200] as a function of fiducial uncertainty. See 
Section VI for details. 



Does this effect future inference? 

Table I summarizes fiducial uncertainty (at 15%), es- 
timation standard deviation, and bias range when using 
4>i^ to estimate the average power in I'^Cf'^ /4, over / e 
[10,200]. Each row corresponds to a different quadratic 
pairing X, F € {T, i?, i?}. The second column shows 
the bias range corresponding to 15% fiducial uncertainty 
taken from figures 4, 5 and 6. We list bias range, ver- 
sus absolute bias, since any baseline bias can be esti- 
mated with simulation and subsequently subtracted from 
any estimate. The third column shows 4(t full sky error 
bars where the estimation standard deviation is extrapo- 
lated from figures 4, 5 and 6 to full sky a by multiplying 
■y/Tsky = 0.126. Notice that, ignoring the bias, the EE 
and EB estimator have the power to constrain the fidu- 
cial uncertainty by a factor of 3. However, accounting for 
bias this constraining power is mitigated, especially for 
estimators EE and TE. In contrast, the EB estimator 



VII. DISCUSSION 

The state-of-the-art estimator of weak lensing, the 
quadratic estimator developed by Hu and Okomoto [3, 4], 
works in part through a delicate cancelation of terms in 
a Taylor expansion of the lensing effect on the CMB. In 
this paper we present two simulation based approaches 
for exploring the nature of this cancelation for both the 
CMB intensity and the polarization fields. In particular, 
we use these two simulation algorithms to analyze the so 
called Ni^'' and A^j^^^ bias for two modifications of the full 
set of quadratic temperature/polarization lensing estima- 
tors: one which incorporates Icnsed rather than unlcnscd 
spectra into the estimator weights to mitigate the effect 
of higher order terms; and one which uses the observed 
lensed CMB fields to correct for the, so called, Ni'^^ bias. 
Our simulation algorithms, which can simulate all higher 

order bias terms N^^'' for j > 0, utilize an extension of the 
FFT techniques developed in [3] . These FFT characteri- 
zations are key to making estimates of the N^^^ fast and 
additionally provide fast (non-stochastic) algorithms for 
approximating the A^/°^ bias using the observed lensed 
CMB fields. 

In Section VI we use our algorithm to analyze the mod- 
ified quadratic temperature/polarization lensing estima- 
tors for future ACTpol/SPTpol experiments. We find 
that the modified estimates do reduce low / estimation 
bias. However this is accomplished by effectively increas- 
ing the magnitude of Ni^^ and Ni"^"^ to the point of can- 
celation when the correct model for Cf"^ is used to gener- 
ate the lensed weights. We also demonstrate, through an 
analysis of estimator bias versus fiducial uncertainty, that 
this cancelation can be sensitive to the fiducial model for 
Cf'^ depending on which estimator one uses: TE and EE 
are most sensitive, EB is least sensitive. For low I estima- 
tion in future ACTpol/SPTpol experiments we conclude 
that the bias in the EB estimator can be effectively ig- 
nored. For the TE and the EE estimators, however, the 
bias does contribute significantly to projected error bars 
and may need to be corrected to give the estimator in- 
ferential power beyond a nominal fiducial uncertainty. 
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Appendix A: FFT algorithms 

In this section we derive the FFT algorithms which al- 
low fast simulation of the fields and ^^^^ as 
described in Section V A. We begin with some notational 
conventions which greatly simplify the subsequent formu- 
las. First let Lpk denote the phase angle of frequency k 
and = ^k+i — fk- Also let la denote the a"^ coordi- 
nate of I. For any field X G {T, E, B} we let 



Xk 



Furthermore, we will utilize a super-script/sub-script no- 
tation to denote multiplication/division by particular 
power spectra. An example serves to illustrate the nota- 
tion: 



1 



= x. 



rC,exp 



\^Tr,EE ) — 



f^TE 



^TT ^EE 
fc,cxp fcjgxp 
^TE^TT 



CI 



(J EE 



Notice that the above denominators always use lensed 
spectra with experimental noise, whereas the numerators 
always use unlensed spectra. In doing so, the formulas 
found in claims 1 through 5 below can be used for fast 
algorithms for the quadratic estimate (j)f^ ■ To obtain 
the corresponding formulas for the modified quadratic 
estimate one simply needs to replace the unlensed 
spectra in the numerator of the above notation, with the 
lensed spectra (but without experimental noise). 

Claim 1 (TT estimator) If X^ and are complex 
functions such that X_k = and Y-k = Y^ then 



d^k 



91% ^fe+t^^fe = -«^a 



d^X 



-ilx 



2lT 

V-X^{x)][Y^{x)] 
+ [V^Y^{x)][X^{x)]); 



d^k 



9i% \ ^k+lYk 



IJb 

2ir 



I 



d^x 



-ilx 



2[y^x:^^^{x)][vVM] 



[X^^MW'^'Y^^-T^ix) 
[Y^,^{xm^V'X^;^{x)] 



where gfl ^ 27r/-[C-,.„ C^^J-i and f^ ^ ^[l 
{k + l)Clli-l-kCl-]. 



Claim 2 (TE estimator) // X^ and Yk are complex 
functions such that X-k = X^ and Y-k = Y^ then 



^9i,kXk+iYk--da 



X [ReiiV^X^ixW^^ix)]*} 
+ [V^Y^^{x)][X^{x)]); 



d^x 



-ilx 



— \g^,^\Xk+iY^--—j 

X {[Y„{xWV''Xl^^;^{x)] 
+ i?e{[y,,,,,(x)] [V^V^X™;- {x)]*} 

+ 2[x^,^(x)][v«v^y™'™(x)] 

+ i?e{4[V'^l™^(x)][V''fr,..(^)]*}); 



27r 



/d^k _ Iglb f 

j^9i,k9i,-k-i^k+iyk --^ ] 

X (iie{2[V"V*X-;S(x)][y^,,,(a:)]*} 

+ Re{2[V'^V'Y:^:S{x)][X^^,,{x)Y} 
+ 3[V«X™,Ja:)][V^y-%(x)] 

+ i?e{[V^y-%(a;)][V«X™,,(x)]*}), 

wheregj% ^ 27r/™[C- ^^^^ C-J"! and/- ^ ^[l-{k+ 
l)cos{2A<f)CZli-l-kC^]. 

Claim 3 (TB estimator) // Xk and Yk are complex 
functions such that X-k = and Y-k = Y^ then 

'^^n™ Y V*— ,7 /" „-il-x 

-^9i,kXk+iYk--ila j 

X Im{[W''X^{x)][%^{x)]*}; 



d^x 



-il-x 



^ \^TB \2\r T^* _ f 

j^\9l,k\ Xk+iY^--—J 2^-. 
X ([YBB,BB(x)][V°V''X™;™(a;)] 

- Re{[Y^s,BB{x)W"^'^^;^ix)]*}) ; 



/ 



^ J ^ 



d'^k II, 
^2-,^y9i,k9i-k-l^k+iyk -~47r 

X ([V»X-,,(a;)][V''y-%(a;)] 



—ila 

^ J 



- Re{[V'^X^^^Ax)][V''Y^%^{x)]*}) , 

wheregf-^ ^ 27r/™[C-,,^,^C-J-i and f[- ^ ^[l- 
i)C™Jsin(2Ay,). 



[l-{k+ 
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Claim 4 (EB estimator) // Xf, and are complex 
functions such that X-k = X^^ and Y-k = Yj^ then 

x7m{[V°lS(a;)][Ye4xr}; 




X ([F,.,BB(a;)][V'^V^X||;||(a;)] 
- Re{[Y,^,s4x)][V-V'XS:Si^)]*}); 

f „EB „EB y V*— f ^'"^ ^-il-X 

J J^9i,k9i,-k-iXk+iY^--— J —e 

- Re{[V^XS,sA=^)WY^E,ss{x)n), 

where gf% = 27r/^f [C,-,,^, C,- J-^ and ff;Z ^ • 
(fe + 0C^fysin(2A^). 

Claim 5 (EE estimator) If Xk and Y^ are complex 



functions such that X-k = and Y-k = Y^ then 

X {Re{[V'^XZ{xW,,{x)Y] 
+ Re{[V'^YSS{x)][X,^{x)]*}y, 




+ [X,,^,,{x)][V-V'Y^^;^^{x)] 

+ Re{[X,,^,,{x)][V'^V'Y^^;^^{x)r} 

+ 2[V^X-,l,,{x)][V'^Y^l,,{x)] 

+ Re{2[V'^K%,{x)][V'YSl^^{x)Y}); 

where gf% ^ 2^//^f [Cf-,.,^,^ CflJ"! and /^f ^ ^[l ■ 
(fe + Z)Cf-,-Z-fcCr]cos(2Av.). 
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